home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / dc-0_2.lha / dc-0.2 / decimal.c < prev    next >
C/C++ Source or Header  |  1993-05-21  |  30KB  |  1,236 lines
  1. /* 
  2.  * Arbitrary precision decimal arithmetic.
  3.  *
  4.  * Copyright (C) 1984 Free Software Foundation, Inc.
  5.  *
  6.  * This program is free software; you can redistribute it and/or modify
  7.  * it under the terms of the GNU General Public License as published by
  8.  * the Free Software Foundation; either version 2, or (at your option)
  9.  * any later version.
  10.  *
  11.  * This program is distributed in the hope that it will be useful,
  12.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.  * GNU General Public License for more details.
  15.  *
  16.  * You should have received a copy of the GNU General Public License
  17.  * along with this program; if not, you can either send email to this
  18.  * program's author (see below) or write to: The Free Software Foundation,
  19.  * Inc.; 675 Mass Ave. Cambridge, MA 02139, USA.
  20.  */
  21.  
  22. /* Some known problems:
  23.  
  24.     Another problem with decimal_div is found when you try to
  25.     divide a number with > scale fraction digits by 1.  The
  26.     expected result is simply truncation, but all sorts of things
  27.     happen instead.  An example is that the result of .99999998/1
  28.     with scale set to 6 is .000001
  29.         
  30.     There are some problems in the behavior of the decimal package
  31.     related to printing and parsing.  The
  32.     printer is weird about very large output radices, tending to want
  33.     to output single ASCII characters for any and all digits (even
  34.     in radices > 127).  The UNIX bc approach is to print digit groups
  35.     separated by spaces.  There is a rather overwrought workaround in
  36.     the function decputc() in bcmisc.c, but it would be better if
  37.     decimal.c got a fix for this.  */
  38.  
  39. /* For stand-alone testing, compile with -DTEST.
  40.    This DTESTable feature defines a `main' function
  41.    which is a simple loop that accepts input of the form
  42.    number space op space number newline
  43.    where op is +, -, *, /, %, p or r,
  44.    and performs the operation and prints the operands and result.
  45.    `p' means print the first number in the radix spec'd by the second.
  46.    `r' means read the first one in the radix specified by the second
  47.    (and print the result in decimal).
  48.    Divide in this test keeps three fraction digits. */
  49.  
  50. #include "decimal.h"
  51.  
  52. #define MAX(a, b) (((a) > (b) ? (a) : (b)))
  53.  
  54. /* Some constant decimal numbers */
  55.  
  56. struct decimal decimal_zero = {0, 0, 0, 0, 0};
  57.  
  58. struct decimal decimal_one = {0, 0, 1, 0, 1};
  59.  
  60. /*** Assumes RADIX is even ***/
  61. struct decimal decimal_half = {0, 1, 0, 0, RADIX / 2};
  62.  
  63. decimal static decimal_add1 (), decimal_sub1 ();
  64. static void add_scaled ();
  65. static int subtract_scaled ();
  66.  
  67. /* Create and return a decimal number that has `before' digits before
  68.    the decimal point and `after' digits after.  The digits themselves are
  69.    initialized to zero.  */
  70.  
  71. decimal
  72. make_decimal (before, after)
  73.      int before, after;
  74. {
  75.   decimal result;
  76.   if (before >= 1<<16)
  77.     {
  78.       decimal_error ("%d too many decimal digits", before);
  79.       return 0;
  80.     }
  81.   if (after >= 1<<15)
  82.     {
  83.       decimal_error ("%d too many decimal digits", after);
  84.       return 0;
  85.     }
  86.   result = (decimal) malloc (sizeof (struct decimal) + before + after - 1);
  87.   result->sign = 0;
  88.   result->before = before;
  89.   result->after = after;
  90.   result->refcnt = 0;
  91.   bzero (result->contents, before + after);
  92.   return result;
  93. }
  94.  
  95. /* Create a copy of the decimal number `b' and return it.  */
  96.  
  97. decimal
  98. decimal_copy (b)
  99.      decimal b;
  100. {
  101.   decimal result = make_decimal (b->before, b->after);
  102.   bcopy (b->contents, result->contents, LENGTH(b));
  103.   result->sign = b->sign;
  104.   return result;
  105. }
  106.  
  107. /* Copy a decimal number `b' but extend or truncate to exactly
  108.    `digits' fraction digits. */
  109.  
  110. static decimal
  111. decimal_copy_1 (b, digits)
  112.      decimal b;
  113.      int digits;
  114. {
  115.   if (digits > b->after)
  116.     {
  117.       decimal result = make_decimal (b->before, digits);
  118.       bcopy (b->contents, result->contents + (digits - (int) b->after), LENGTH(b));
  119.       return result;
  120.     }
  121.   else
  122.     return decimal_round_digits (b, digits);
  123. }
  124.  
  125. /* flush specified number `digits' of trailing fraction digits,
  126.    and flush any trailing fraction zero digits exposed after they are gone.
  127.    The number `b' is actually modified; no new storage is allocated.
  128.    That is why this is not global.  */
  129.  
  130. static void
  131. flush_trailing_digits (b, digits)
  132.      decimal b;
  133.      int digits;
  134. {
  135.   int flush = digits;
  136.   int maxdig = b->after;
  137.  
  138.   while (flush < maxdig && !b->contents [flush])
  139.     flush++;
  140.  
  141.   if (flush)
  142.     {
  143.       int i;
  144.  
  145.       b->after -= flush;
  146.       for (i = 0; i < LENGTH (b); i++)
  147.     b->contents[i] = b->contents[flush + i];
  148.     }
  149.  
  150. }
  151.  
  152. /* Return nonzero integer if the value of decimal number `b' is zero.  */
  153.  
  154. int
  155. decimal_zerop (b)
  156.      decimal b;
  157. {
  158.   return !LENGTH(b);
  159. }
  160.  
  161. /* Compare two decimal numbers arithmetically.
  162.    The value is < 0 if b1 < b2, > 0 if b1 > b2, 0 if b1 = b2.
  163.    This is the same way that `strcmp' reports the result of comparing
  164.    strings.  */ 
  165.  
  166. int
  167. decimal_compare (b1, b2)
  168.      decimal b1, b2;
  169. {
  170.   int l1, l2;
  171.   char *p1, *p2, *s1, *s2;
  172.   int i;
  173.  
  174.   /* If signs differ, deduce result from the signs */
  175.  
  176.   if (b2->sign && !b1->sign) return 1;
  177.   if (b1->sign && !b2->sign) return -1;
  178.  
  179.   /* If same sign but number of nonfraction digits differs,
  180.      the one with more of them is farther from zero.  */
  181.  
  182.   if (b1->before != b2->before)
  183.     if (b1->sign)
  184.       return (int) (b2->before - b1->before);
  185.     else
  186.       return (int) (b1->before - b2->before);
  187.  
  188.   /* Else compare the numbers digit by digit from high end */
  189.   l1 = LENGTH(b1);
  190.   l2 = LENGTH(b2);  
  191.   s1 = b1->contents;        /* Start of number -- don't back up digit pointer past here */
  192.   s2 = b2->contents;
  193.   p1 = b1->contents + l1;    /* Scanning pointer, for fetching digits.  */
  194.   p2 = b2->contents + l2;
  195.   for (i = MAX(l1, l2); i >= 0; i--)
  196.     {
  197.       int r = ((p1 != s1) ? *--p1 : 0) - ((p2 != s2) ? *--p2 : 0);
  198.       if (r)
  199.     return b1->sign ? -r : r;
  200.     }
  201.   return 0;
  202. }
  203.  
  204. /* Return the number of digits stored in decimal number `b' */
  205.  
  206. int
  207. decimal_length (b)
  208.      decimal b;
  209. {
  210.   return LENGTH(b);
  211. }
  212.  
  213. /* Return the number of fraction digits stored in decimal number `b'.  */
  214.  
  215. int
  216. decimal_after (b)
  217.      decimal b;
  218. {
  219.   return b->after;
  220. }
  221.  
  222. /* Round decimal number `b' to have only `digits' fraction digits.
  223.    Result is rounded to nearest unit in the last remaining digit.
  224.    Return the result, another decimal number.  */
  225.  
  226. decimal
  227. decimal_round_digits (b, digits)
  228.      decimal b;
  229.      int digits;
  230. {
  231.   decimal result;
  232.   int old;
  233.  
  234.   if (b->after <= digits) return decimal_copy (b);
  235.  
  236.   if (digits < 0)
  237.     {
  238.       decimal_error ("request to keep negative number of digits %d", digits);
  239.       return decimal_copy (b);
  240.     }
  241.  
  242.   result = make_decimal (b->before + 1, b->after);
  243.   result->sign = b->sign;
  244.   bcopy (b->contents, result->contents, LENGTH(b));
  245.  
  246.   old = result->after;
  247.  
  248.   /* Add .5 * last place to keep, so that we round rather than truncate */
  249.   /* Note this ignores sign of result, so if result is negative
  250.      it is subtracting */
  251.  
  252.   add_scaled (result, DECIMAL_HALF, 1, old - digits - 1);
  253.  
  254.   /* Flush desired digits, and any trailing zeros exposed by them.  */
  255.  
  256.   flush_trailing_digits (result, old - digits);
  257.  
  258.   /* Flush leading digits -- always is one, unless was a carry into it */
  259.  
  260.   while (result->before > 0
  261.      && result->contents[LENGTH(result) - 1] == 0)
  262.     result->before--;
  263.  
  264.   return result;
  265. }
  266.  
  267. /* Truncate decimal number `b' to have only `digits' fraction digits.
  268.    Any fraction digits in `b' beyond that are dropped and ignored.
  269.    Truncation is toward zero.
  270.    Return the result, another decimal number.  */
  271.  
  272. decimal
  273. decimal_trunc_digits (b, digits)
  274.      decimal b;
  275.      int digits;
  276. {
  277.   decimal result = decimal_copy (b);
  278.   int old = result->after;
  279.  
  280.   if (old <= digits) return result;
  281.  
  282.   if (digits < 0)
  283.     {
  284.       decimal_error ("request to keep negative number of digits %d", digits);
  285.       return result;
  286.     }
  287.  
  288.   flush_trailing_digits (result, old - digits);
  289.  
  290.   return result;
  291. }
  292.  
  293. /* Return the fractional part of decimal number `b':
  294.    that is, `b' - decimal_trunc_digits (`b') */
  295.  
  296. decimal
  297. decimal_fraction (b)
  298.      decimal b;
  299. {
  300.   decimal result = make_decimal (0, b->after);
  301.   bcopy (b->contents, result->contents, b->after);
  302.   return result;
  303. }
  304.  
  305. /* return an integer whose value is that of decimal `b', sans its fraction.  */
  306.  
  307. int
  308. decimal_to_int (b)
  309.      decimal b;
  310. {
  311.   int result = 0;
  312.   int i;
  313.   int end = b->after;
  314.  
  315.   for (i = LENGTH(b) - 1; i >= end; i--)
  316.     {
  317.       result *= RADIX;
  318.       result += b->contents[i];
  319.     }
  320.   return result;
  321. }
  322.  
  323. /* return a decimal whose value is the integer i.  */
  324.  
  325. decimal
  326. decimal_from_int (i)
  327.      int i;
  328. {
  329.   int log, tem;
  330.   decimal result;
  331.  
  332.   for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX);
  333.  
  334.   result = make_decimal (log, 0);
  335.  
  336.   for (log = 0, tem = (i > 0 ? i : - i); tem; log++, tem /= RADIX)
  337.     result->contents[log] = tem % RADIX;
  338.  
  339.   if (i < 0) result->sign = 1;
  340.   return result;
  341. }
  342.  
  343. /* Return (as an integer) the result of dividing decimal number `b' by
  344.    integer `divisor'. 
  345.    This is used in printing decimal numbers in other radices. */
  346.  
  347. int
  348. decimal_int_rem (b, divisor)
  349.      decimal b;
  350.      int divisor;
  351. {
  352.   int len = LENGTH(b);
  353.   int end = b->after;
  354.   int accum = 0;
  355.   int i;
  356.  
  357.   for (i = len - 1; i >= end; i--)
  358.     {
  359.       accum %= divisor;
  360.       accum *= RADIX;
  361.       accum += b->contents[i];
  362.     }
  363.   return accum % divisor;
  364. }
  365.  
  366. /* Convert digit `digit' to a character and output it by calling
  367.    `charout' with it as arg. */
  368.  
  369. static void
  370. print_digit (digit, charout)
  371.      int digit;
  372.      void (*charout) ();
  373. {
  374.   if (digit < 10)
  375.     charout ('0' + digit);
  376.   else
  377.     charout ('A' + digit - 10);
  378. }
  379.  
  380. /* print decimal number `b' in radix `radix', assuming it is an integer.
  381.    `r' is `radix' expressed as a decimal number. */
  382.  
  383. static
  384. decimal_print_1 (b, r, radix, charout)
  385.      decimal b, r;
  386.      int radix;
  387.      void (*charout) ();
  388. {
  389.   int digit = decimal_int_rem (b, radix);
  390.   decimal rest = decimal_div (b, r, 0);
  391.  
  392.   if (!decimal_zerop (rest))
  393.     decimal_print_1 (rest, r, radix, charout);
  394.  
  395.   print_digit (digit, charout);
  396.  
  397.   free (rest);
  398. }
  399.  
  400. /* User entry: print decimal number `b' in radix `radix' (an integer),
  401.    outputting characters by calling `charout'.  */
  402.  
  403. void
  404. decimal_print (b, charout, radix)
  405.      decimal b;
  406.      void (*charout) ();
  407.      int radix;
  408. {
  409.   if (b->sign) charout ('-');
  410.  
  411.   if (radix == RADIX)
  412.     {
  413.       /* decimal output => just print the digits, inserting a point in
  414.      the proper place.  */ 
  415.       int i;
  416.       int before = b->before;
  417.       int len = before + b->after;
  418.       for (i = 0; i < len; i++)
  419.     {
  420.       if (i == before) charout ('.');
  421.       /* Broken if RADIX /= 10
  422.          charout ('0' + b->contents [len - 1 - i]); */
  423.       print_digit (b->contents [len - 1 - i], charout);
  424.     }
  425.       if (!len)
  426.     charout ('0');
  427.     }
  428.   else
  429.     {
  430.       /* nonstandard radix: must use multiply and divide to determine the
  431.      digits of the number in that radix.  */
  432.  
  433.       int i;
  434.       extern double log10 ();
  435.       /* Compute the number of fraction digits we want to have in the
  436.          new radix.  They should contain the same amount of
  437.          information as the decimal digits we have.  */
  438.       int nfrac = (b->after / log10 ((double) radix) + .99);
  439.       decimal r = decimal_from_int (radix);
  440.       decimal intpart = decimal_trunc_digits (b, 0);
  441.  
  442.       /* print integer part */
  443.       decimal_print_1 (intpart, r, radix, charout);
  444.       free (intpart);
  445.  
  446.       /* print fraction part */
  447.       if (nfrac)
  448.     {
  449.           decimal tem1, tem2;
  450.       tem1 = decimal_fraction (b);
  451.       charout ('.');
  452.       /* repeatedly multiply by `radix', print integer part as one digit,
  453.          and flush the integer part.  */
  454.       for (i = 0; i < nfrac; i++)
  455.         {
  456.           tem2 = decimal_mul (tem1, r);
  457.           free (tem1);
  458.           print_digit (decimal_to_int (tem2), charout);
  459.           tem1 = decimal_fraction (tem2);
  460.           free (tem2);
  461.         }
  462.       free (tem1);
  463.     }
  464.       free (r);
  465.     }
  466. }
  467.  
  468. static int
  469. decode_digit (digitchar)
  470.      char digitchar;
  471. {
  472.   if ('0' <= digitchar && digitchar <= '9')
  473.     return digitchar - '0';
  474.   if ('a' <= digitchar && digitchar <= 'z')
  475.     return digitchar - 'a' + 10;
  476.   if ('A' <= digitchar && digitchar <= 'Z')
  477.     return digitchar - 'A' + 10;
  478.   return -1;
  479. }
  480.  
  481. /* Parse string `s' into a number using radix `radix'
  482.    and return result as a decimal number.  */
  483.  
  484. decimal
  485. decimal_parse (s, radix)
  486.      char *s;
  487.      int radix;
  488. {
  489.   int i, len, before = -1;
  490.   char *p;
  491.   char c;
  492.   decimal result;
  493.   int negative = 0;
  494.   int excess_digit = 0;
  495.  
  496.   if (*s == '-')
  497.     {
  498.       s++;
  499.       negative = 1;
  500.     }
  501.  
  502.   /* First scan for valid characters.
  503.      Count total num digits, and count num before the decimal point.  */
  504.  
  505.   p = s;
  506.   i = 0;
  507.   while (c = *p++)
  508.     {
  509.       if (c == '.')
  510.         {
  511.       if (before >= 0)
  512.         decimal_error ("two decimal points in %s", s);
  513.           before = i;
  514.     }
  515.       else if (c == '0' && !i && before < 0)
  516.     s++;   /* Discard leading zeros */
  517.       else if (decode_digit (c) >= 0)
  518.     {
  519.       i++;
  520.       if (decode_digit (c) > RADIX)
  521.         excess_digit = 1;
  522.     }
  523.       else
  524.     decimal_error ("invalid number %s", s);
  525.     }
  526.  
  527.   len = i;
  528.   if (before < 0) before = i;
  529.  
  530.   p = s;
  531.  
  532.   /* Now parse those digits */
  533.  
  534.   if (radix != RADIX || excess_digit)
  535.     {
  536.       decimal r = decimal_from_int (radix);
  537.       extern double log10 ();
  538.       int digits = (len - before) * log10 ((double) radix) + .99;
  539.       result = decimal_copy (DECIMAL_ZERO);
  540.  
  541.       /* Parse all the digits into an integer, ignoring decimal point,
  542.      by multiplying by `radix'.  */
  543.  
  544.       while (i > 0 && (c = *p++))
  545.     {
  546.       if (c != '.')
  547.         {
  548.           decimal newdig = decimal_from_int (decode_digit (c));
  549.           decimal prod = decimal_mul (result, r);
  550.           decimal newresult = decimal_add (newdig, prod);
  551.  
  552.           free (newdig);  free (prod);  free (result);
  553.           result = newresult;
  554.           i--;
  555.         }
  556.     }
  557.  
  558.       /* Now put decimal point in right place
  559.      by dividing by `radix' once for each digit
  560.      that really should have followed the decimal point.  */
  561.  
  562.       for (i = before; i < len; i++)
  563.     {
  564.       decimal newresult = decimal_div (result, r, digits);
  565.       free (result);
  566.       result = newresult;
  567.     }
  568.       free (r);
  569.     }
  570.   else
  571.     {
  572.       /* radix is standard - just copy the digits into a decimal number.  */
  573.  
  574.       int tem;
  575.       result = make_decimal (before, len - before);
  576.  
  577.       while (i > 0 && (c = *p++))
  578.     {
  579.       if ((c != '.') &&
  580.           ((tem = decode_digit (c)) >= 0))
  581.         result->contents [--i] = tem;
  582.     }
  583.     }
  584.  
  585.   if (negative) result->sign = 1;
  586.   flush_trailing_digits (result, 0);
  587.   return result;
  588. }
  589.  
  590. /* Add b1 and b2, considering their signs */
  591.  
  592. decimal
  593. decimal_add (b1, b2)
  594.      decimal b1, b2;
  595. {
  596.   decimal v;
  597.  
  598.   if (b1->sign != b2->sign)
  599.     v = decimal_sub1 (b1, b2);
  600.   else
  601.     v = decimal_add1 (b1, b2);
  602.   if (b1->sign && !decimal_zerop (v))
  603.     v->sign = !v->sign;
  604.   return v;
  605. }
  606.  
  607. /* Add b1 and minus b2, considering their signs */
  608.  
  609. decimal
  610. decimal_sub (b1, b2)
  611.      decimal b1, b2;
  612. {
  613.   decimal v;
  614.  
  615.   if (b1->sign != b2->sign)
  616.     v = decimal_add1 (b1, b2);
  617.   else
  618.     v = decimal_sub1 (b1, b2);
  619.   if (b1->sign && !decimal_zerop (v))
  620.     v->sign = !v->sign;
  621.   return v;
  622. }
  623.  
  624. /* Return the negation of b2.  */
  625.  
  626. decimal
  627. decimal_neg (b2)
  628.      decimal b2;
  629. {
  630.   decimal v = decimal_copy (b2);
  631.  
  632.   if (!decimal_zerop (v))
  633.     v->sign = !v->sign;
  634.   return v;
  635. }
  636.  
  637. /* add magnitudes of b1 and b2, ignoring their signs. */
  638.  
  639. static decimal
  640. decimal_add1 (b1, b2)
  641.      decimal b1, b2;
  642. {
  643.   int before = MAX (b1->before, b2->before);
  644.   int after = MAX (b1->after, b2->after);
  645.  
  646.   int len = before+after+1;
  647.   decimal result = make_decimal (before+1, after);
  648.  
  649.   int i;
  650.   char *s1 = b1->contents;
  651.   char *s2 = b2->contents;
  652.   char *p1 = s1 + b1->after - after;
  653.   char *p2 = s2 + b2->after - after;
  654.   char *e1 = s1 + b1->before + b1->after;
  655.   char *e2 = s2 + b2->before + b2->after;
  656.   char *pr = result->contents;
  657.   int accum = 0;
  658.  
  659.   for (i = 0; i < len; i++, p1++, p2++)
  660.     {
  661.       accum /= RADIX;
  662.       if (p1 >= s1 && p1 < e1) accum += *p1;
  663.       if (p2 >= s2 && p2 < e2) accum += *p2;
  664.       *pr++ = accum % RADIX;
  665.     }
  666.   if (!accum)
  667.     (result->before)--;
  668.  
  669.   flush_trailing_digits (result, 0);
  670.  
  671.   return result;
  672. }
  673.  
  674. /* subtract magnitude of b2 from that or b1, returning signed decimal
  675.    number. */ 
  676.  
  677. static decimal
  678. decimal_sub1 (b1, b2)
  679.      decimal b1, b2;
  680. {
  681.   int before = MAX (b1->before, b2->before);
  682.   int after = MAX (b1->after, b2->after);
  683.  
  684.   int len = before+after;
  685.   decimal result = make_decimal (before, after);
  686.  
  687.   int i;
  688.   char *s1 = b1->contents;
  689.   char *s2 = b2->contents;
  690.   char *p1 = s1 + b1->after - after;
  691.   char *p2 = s2 + b2->after - after;
  692.   char *e1 = s1 + b1->before + b1->after;
  693.   char *e2 = s2 + b2->before + b2->after;
  694.   char *pr = result->contents;
  695.   int accum = 0;
  696.  
  697.   for (i = 0; i < len; i++, p1++, p2++)
  698.     {
  699.       if (p1 >= s1 && p1 < e1) accum += *p1;
  700.       if (p2 >= s2 && p2 < e2) accum -= *p2;
  701.       if (accum < 0 && accum % RADIX)
  702.         *pr = RADIX - (- accum) % RADIX;
  703.       else
  704.     *pr = accum % RADIX;
  705.       accum -= *pr++;
  706.       accum /= RADIX;
  707.     }
  708.  
  709.   /* If result is negative, subtract it from RADIX**length
  710.      so that we get the right digits for sign-magnitude
  711.      rather than RADIX-complement */
  712.  
  713.   if (accum)
  714.     {
  715.       result->sign = 1;
  716.       pr = result->contents;
  717.       accum = 0;
  718.       for (i = 0; i < len; i++)
  719.     {
  720.       accum -= *pr;
  721.       if (accum)
  722.         *pr = accum + RADIX;
  723.       else
  724.         *pr = 0;
  725.       accum -= *pr++;
  726.       accum /= RADIX;
  727.     }
  728.     }
  729.  
  730.   /* flush leading nonfraction zero digits */
  731.  
  732.   while (result->before && *--pr == 0)
  733.     (result->before)--;
  734.  
  735.   flush_trailing_digits (result, 0);
  736.  
  737.   return result;
  738. }
  739.  
  740. /* multiply b1 and b2 keeping `digits' fraction digits */
  741.  
  742. decimal
  743. decimal_mul_rounded (b1, b2, digits)
  744.      decimal b1, b2;
  745.      int digits;
  746. {
  747.   decimal tem = decimal_mul (b1, b2);
  748.   decimal result = decimal_round_digits (tem, digits);
  749.   free (tem);
  750.   return result;
  751. }
  752.  
  753. /* multiply b1 and b2 keeping the right number of fraction digits
  754.    for the `dc' program with precision = `digits'.  */
  755.  
  756. decimal
  757. decimal_mul_dc (b1, b2, digits)
  758.      decimal b1, b2;
  759.      int digits;
  760. {
  761.   decimal tem = decimal_mul (b1, b2);
  762.   decimal result
  763.     = decimal_round_digits (tem, MAX (digits, MAX (b1->after, b2->after)));
  764.   free (tem);
  765.   return result;
  766. }
  767.  
  768. /* multiply b1 and b2 as decimal error-free values;
  769.    keep LENGTH(b1) plus LENGTH(b2) significant figures. */
  770.  
  771. decimal
  772. decimal_mul (b1, b2)
  773.      decimal b1, b2;
  774. {
  775.   decimal result = make_decimal (b1->before + b2->before, b1->after + b2->after);
  776.   int i;
  777.   int length2 = LENGTH(b2);
  778.   char *pr;
  779.  
  780.   for (i = 0; i < length2; i++)
  781.     add_scaled (result, b1, b2->contents[i], i);
  782.  
  783.   /* flush leading nonfraction zero digits */
  784.  
  785.   pr = result->contents + LENGTH(result);
  786.   while (result->before && *--pr == 0)
  787.     (result->before)--;
  788.  
  789.   flush_trailing_digits (result, 0);   /* flush trailing zeros */
  790.  
  791.   /* Set sign properly */
  792.  
  793.   if (b1->sign != b2->sign && LENGTH(result))
  794.     result->sign = 1;
  795.  
  796.   return result;
  797. }
  798.  
  799. /* Modify decimal number `into' by adding `from',
  800.    multiplied by `factor' (which should be nonnegative and less than RADIX)
  801.    and shifted left `scale' digits at the least significant end. */
  802.  
  803. static void
  804. add_scaled (into, from, factor, scale)
  805.      decimal into, from;
  806.      int factor, scale;
  807. {
  808.   char *pf = from->contents;
  809.   char *pi = into->contents + scale;
  810.   int lengthf = LENGTH(from);
  811.   int lengthi = LENGTH(into) - scale;
  812.   
  813.   int accum = 0;
  814.   int i;
  815.  
  816.   for (i = 0; i < lengthi; i++)
  817.     {
  818.       accum /= RADIX;
  819.       if (i < lengthf)
  820.         accum += *pf++ * factor;
  821.       accum += *pi;
  822.       *pi++ = accum % RADIX;
  823.     }
  824. }
  825.  
  826. /* Divide decimal number `b1' by `b2', keeping at most `digits'
  827.    fraction digits. 
  828.    Returns the result as a decimal number.
  829.  
  830.    When division is not exact, the quotient is truncated toward zero.  */
  831.  
  832. decimal
  833. decimal_div (b1, b2, digits)
  834.      decimal b1, b2;
  835.      int digits;
  836. {
  837.   decimal result = make_decimal (MAX(1, (int) (1 + b1->before - b2->before)), digits);
  838.  
  839.   /* b1copy holds what is left of the dividend,
  840.      that is not accounted for by the quotient digits already known */
  841.  
  842.   decimal b1copy = decimal_copy_1 (b1, b2->after + digits);
  843.   int length1 = LENGTH(b1copy);
  844.   int length2 = LENGTH(b2);
  845.   int lengthr = LENGTH(result);
  846.   int i;
  847.  
  848.   /* leading_divisor_digits contains the first two divisor digits, as
  849.      an integer */ 
  850.  
  851.   int leading_divisor_digits = b2->contents[length2-1]*RADIX;
  852.   if (length2 > 1)
  853.     leading_divisor_digits += b2->contents[length2-2];
  854.  
  855.   if (decimal_zerop (b2))
  856.     {
  857.       decimal_error ("divisor is zero", 0);
  858.       return decimal_copy (DECIMAL_ZERO);
  859.     }
  860.  
  861.   if (lengthr <= (length1 - length2))
  862.     abort();         /* My reasoning says this cannot happen, I hope */
  863.  
  864.   for (i = length1 - length2; i >= 0; i--)
  865.     {
  866.       /* Guess the next quotient digit (in order of decreasing significance)
  867.      using integer division */
  868.  
  869.       int guess;
  870.       int trial_dividend = b1copy->contents[length2+i-1]*RADIX;
  871.       if (i != length1 - length2)
  872.     trial_dividend += b1copy->contents[length2+i]*RADIX*RADIX;
  873.       if (length2 + i > 1)
  874.     trial_dividend += b1copy->contents[length2+i-2];
  875.  
  876.       guess = trial_dividend / leading_divisor_digits;
  877.  
  878.       /* Remove the quotient times this digit from the dividend left */
  879.       /* We may find that the quotient digit is too large,
  880.      when we consider the entire divisor.
  881.      Then we decrement the quotient digit and add the divisor back in */
  882.  
  883.       if (guess && 0 > subtract_scaled (b1copy, b2, guess, i))
  884.     {
  885.       guess--;
  886.       add_scaled (b1copy, b2, 1, i);
  887.     }
  888.  
  889.       if (guess >= RADIX)
  890.     {
  891.       result->contents[i + 1] += guess / RADIX;
  892.       guess %= RADIX;
  893.     }
  894.       result->contents[i] = guess;
  895.     }
  896.  
  897.   free (b1copy);
  898.  
  899.   result->sign = (b1->sign != b2->sign);
  900.  
  901.   /* flush leading nonfraction zero digits */
  902.  
  903.   {
  904.     char *pr = result->contents + lengthr;
  905.     while (result->before && *--pr == 0)
  906.       (result->before)--;
  907.   }
  908.  
  909.   flush_trailing_digits (result, 0);    /* Flush trailing zero fraction digits */
  910.  
  911.   return result;
  912. }
  913.  
  914. /* The remainder for the above division.
  915.    Same as `b1' - (`b1' / `b2') * 'b2'.
  916.    Note that the value depends on the number of fraction digits
  917.    that were kept in computing `b1' / `b2';
  918.    the argument `digits' specifies this.
  919.  
  920.    The remainder has the same sign as the dividend.
  921.    The divisor's sign is ignored.  */
  922.  
  923. decimal
  924. decimal_rem (b1, b2, digits)
  925.      decimal b1, b2;
  926.      int digits;
  927. {
  928.   decimal b1copy = decimal_copy_1 (b1, b2->after + digits);
  929.   int length1 = LENGTH(b1copy);
  930.   int length2 = LENGTH(b2);
  931.   int i;
  932.  
  933.   int leading_divisor_digits = b2->contents[length2-1]*RADIX;
  934.  
  935.   if (length2 > 1)
  936.     leading_divisor_digits += b2->contents[length2-2];
  937.  
  938.   if (decimal_zerop (b2))
  939.     {
  940.       decimal_error ("divisor is zero", 0);
  941.       return decimal_copy (DECIMAL_ZERO);
  942.     }
  943.  
  944.   /* Do like division, above, but throw away the quotient.
  945.      Keep only the final `rest of dividend', which becomes the remainder.  */
  946.  
  947.   for (i = length1 - length2; i >= 0; i--)
  948.     {
  949.       int guess;
  950.       int trial_dividend = b1copy->contents[length2+i-1]*RADIX;
  951.       if (i != length1 - length2)
  952.     trial_dividend += b1copy->contents[length2+i]*RADIX*RADIX;
  953.       if (length2 + i > 1)
  954.     trial_dividend += b1copy->contents[length2+i-2];
  955.  
  956.       guess = trial_dividend / leading_divisor_digits;
  957.  
  958.       if (guess && 0 > subtract_scaled (b1copy, b2, guess, i))
  959.     {
  960.       guess--;
  961.       add_scaled (b1copy, b2, 1, i);
  962.     }
  963.       /* No need to check whether guess exceeds RADIX
  964.      since we are not saving guess.  */
  965.     }
  966.  
  967.   /* flush leading nonfraction zero digits */
  968.  
  969.   {
  970.     char *pr = b1copy->contents + length1;
  971.     while (b1copy->before && *--pr == 0)
  972.       (b1copy->before)--;
  973.   }
  974.  
  975.   flush_trailing_digits (b1copy, 0);
  976.   return b1copy;
  977. }
  978.  
  979. /* returns negative number if we chose factor too large */
  980.  
  981. static int
  982. subtract_scaled (into, from, factor, scale)
  983.      decimal into, from;
  984.      int factor, scale;
  985. {
  986.   char *pf = from->contents;
  987.   char *pi = into->contents + scale;
  988.   int lengthf = LENGTH(from);
  989.   int lengthi = LENGTH(into) - scale;
  990.   int accum = 0;
  991.   int i;
  992.  
  993.   for (i = 0; i < lengthi && i <= lengthf; i++)
  994.     {
  995.       if (i < lengthf)
  996.         accum -= *pf++ * factor;
  997.       accum += *pi;
  998.       if (accum < 0 && accum % RADIX)
  999.         *pi = RADIX - (- accum) % RADIX;
  1000.       else
  1001.     *pi = accum % RADIX;
  1002.       accum -= *pi++;
  1003.       accum /= RADIX;
  1004.     }
  1005.   return accum;
  1006. }
  1007.  
  1008. /* Return the square root of decimal number D, using Newton's method.
  1009.    Number of fraction digits returned is max of FRAC_DIGITS
  1010.    and D's number of fraction digits.  */
  1011.  
  1012. decimal
  1013. decimal_sqrt (d, frac_digits)
  1014.      decimal d;
  1015.      int frac_digits;
  1016. {
  1017.   decimal guess;
  1018.   int notdone = 1;
  1019.  
  1020.   if (decimal_zerop (d)) return d;
  1021.   if (d->sign)
  1022.     {
  1023.       decimal_error ("square root argument negative", 0);
  1024.       return decimal_copy (DECIMAL_ZERO);
  1025.     }
  1026.  
  1027.   frac_digits = MAX (frac_digits, d->after);
  1028.  
  1029.   /* Compute an initial guess by taking the square root 
  1030.      of a nearby power of RADIX.  */
  1031.  
  1032.   if (d->before)
  1033.     {
  1034.       guess = make_decimal ((d->before + 1) / 2, 0);
  1035.       guess->contents[guess->before - 1] = 1;
  1036.     }
  1037.   else
  1038.     {
  1039.       /* Arg is less than 1; compute nearest power of RADIX */
  1040.       char *p = d->contents + LENGTH(d);
  1041.       char *sp = p;
  1042.  
  1043.       while (!*--p);    /* Find most significant nonzero digit */
  1044.       if (sp - p == 1)
  1045.     {
  1046.       /* Arg is bigger than 1/RADIX; use 1 as a guess */
  1047.       guess = decimal_copy (DECIMAL_ONE);
  1048.     }
  1049.       else
  1050.     {
  1051.       guess = make_decimal (0, (sp - p) / 2);
  1052.       guess->contents[0] = 1;
  1053.     }
  1054.     }
  1055.  
  1056.   /* Iterate doing guess = (guess + d/guess) / 2  */
  1057.  
  1058.   while (notdone)
  1059.     {
  1060.       decimal tem1 = decimal_div (d, guess, frac_digits + 1);
  1061.       decimal tem2 = decimal_add (guess, tem1);
  1062.       decimal tem3 = decimal_mul_rounded (tem2, DECIMAL_HALF, frac_digits);
  1063.       notdone = decimal_compare (guess, tem3);
  1064.       free (tem1);
  1065.       free (tem2);
  1066.       free (guess);
  1067.       guess = tem3;
  1068.       if (decimal_zerop (guess)) return guess;  /* Avoid divide-by-zero */
  1069.     }
  1070.  
  1071.   return guess;
  1072. }
  1073.  
  1074. /* Raise decimal number `base' to power of integer part of decimal
  1075.    number `expt'.
  1076.    This function depends on using radix 10.
  1077.    It is too hard to write it to work for any value of RADIX,
  1078.    so instead it is simply not available if RADIX is not ten.  */
  1079.  
  1080. #if !(RADIX - 10)
  1081.  
  1082. decimal
  1083. decimal_expt (base, expt, frac_digits)
  1084.      decimal base, expt;
  1085.      int frac_digits;
  1086. {
  1087.   decimal accum = decimal_copy (DECIMAL_ONE);
  1088.   decimal basis1 = base;
  1089.   int digits = expt->before;
  1090.   int dig = 0;                /* Expt digit being processed */
  1091.  
  1092.   if (expt->sign)
  1093.   /* If negative power, take reciprocal first thing
  1094.      so that fraction digit truncation won't destroy
  1095.      what will ultimately be nonfraction digits.  */
  1096.     basis1 = decimal_div (DECIMAL_ONE, base, frac_digits);
  1097.   while (dig < digits)
  1098.     {
  1099.       decimal basis2, basis4, basis8, basis10;
  1100.       int thisdigit = expt->contents[expt->after + dig];
  1101.  
  1102.       /* Compute factors to multiply in for each bit of this digit */
  1103.  
  1104.       basis2 = decimal_mul_rounded (basis1, basis1, frac_digits);
  1105.       basis4 = decimal_mul_rounded (basis2, basis2, frac_digits);
  1106.       basis8 = decimal_mul_rounded (basis4, basis4, frac_digits);
  1107.  
  1108.       /* Now accumulate the factors this digit value selects */
  1109.  
  1110.       if (thisdigit & 1)
  1111.     {
  1112.       decimal accum1 = decimal_mul_rounded (accum, basis1, frac_digits);
  1113.       free (accum);
  1114.       accum = accum1;
  1115.     }
  1116.  
  1117.       if (thisdigit & 2)
  1118.     {
  1119.       decimal accum1 = decimal_mul_rounded (accum, basis2, frac_digits);
  1120.       free (accum);
  1121.       accum = accum1;
  1122.     }
  1123.  
  1124.       if (thisdigit & 4)
  1125.     {
  1126.       decimal accum1 = decimal_mul_rounded (accum, basis4, frac_digits);
  1127.       free (accum);
  1128.       accum = accum1;
  1129.     }
  1130.  
  1131.       if (thisdigit & 8)
  1132.     {
  1133.       decimal accum1 = decimal_mul_rounded (accum, basis8, frac_digits);
  1134.       free (accum);
  1135.       accum = accum1;
  1136.     }
  1137.  
  1138.       /* If there are further digits, compute the basis1 for the next digit */
  1139.  
  1140.       if (++dig < digits)
  1141.     basis10 = decimal_mul_rounded (basis2, basis8, frac_digits);
  1142.  
  1143.       /* Free intermediate results */
  1144.  
  1145.       if (basis1 != base) free (basis1);
  1146.       free (basis2);
  1147.       free (basis4);
  1148.       free (basis8);
  1149.       basis1 = basis10;
  1150.     }
  1151.   return accum;
  1152. }
  1153. #endif
  1154.  
  1155. #ifdef TEST
  1156.  
  1157. fputchar (c)
  1158.      char c;
  1159. {
  1160.   putchar (c);
  1161. }
  1162.  
  1163. /* Top level that can be used to test the arithmetic functions */
  1164.  
  1165. main ()
  1166. {
  1167.   char s1[40], s2[40];
  1168.   decimal b1, b2, b3;
  1169.   char c;
  1170.  
  1171.   while (1)
  1172.     {
  1173.       scanf ("%s %c %s", s1, &c, s2);
  1174.       b1 = decimal_parse (s1, RADIX);
  1175.       b2 = decimal_parse (s2, RADIX);
  1176.       switch (c)
  1177.     {
  1178.     default:
  1179.       c = '+';
  1180.     case '+':
  1181.       b3 = decimal_add (b1, b2);
  1182.       break;
  1183.     case '*':
  1184.       b3 = decimal_mul (b1, b2);
  1185.       break;
  1186.         case '/':
  1187.       b3 = decimal_div (b1, b2, 3);
  1188.       break;
  1189.     case '%':
  1190.       b3 = decimal_rem (b1, b2, 3);
  1191.       break;
  1192.         case 'p':
  1193.       decimal_print (b1, fputchar, RADIX);
  1194.       printf (" printed in base %d is ", decimal_to_int (b2));
  1195.       decimal_print (b1, fputchar, decimal_to_int (b2));
  1196.       printf ("\n");
  1197.       continue;
  1198.     case 'r':
  1199.       printf ("%s read in base %d is ", s1, decimal_to_int (b2));
  1200.       decimal_print (decimal_parse (s1, decimal_to_int (b2)), fputchar, RADIX);
  1201.       printf ("\n");
  1202.       continue;
  1203.     }
  1204.       decimal_print (b1, fputchar, RADIX);
  1205.       printf (" %c ", c);
  1206.       decimal_print (b2, fputchar, RADIX);
  1207.       printf (" = ");
  1208.       decimal_print (b3, fputchar, RADIX);
  1209.       printf ("\n");
  1210.     }
  1211. }
  1212.  
  1213. decimal_error (s1, s2)
  1214.      char *s1, *s2;
  1215. {
  1216.   printf ("\n");
  1217.   printf (s1, s2);
  1218.   printf ("\n");
  1219. }
  1220.  
  1221. static void
  1222. pbi (b)
  1223.       int b;
  1224. {
  1225.   decimal_print ((decimal) b, fputchar, RADIX);
  1226. }
  1227.  
  1228. static void
  1229. pb (b)
  1230.       decimal b;
  1231. {
  1232.   decimal_print (b, fputchar, RADIX);
  1233. }
  1234.  
  1235. #endif
  1236.